AIAA-2001-0986 



Computation of an Underexpanded 3-D 
Rectangular Jet by the CE / SE Method 

Ching Y. Loh, Ananda Himansu, and Xiao Y. Wang 
Taitech Inc., Cleveland, Ohio 

Philip C.E. Jorgenson 

Glenn Research Center, Cleveland, Ohio 


December 2000 


The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated to 
the advancement of aeronautics and space 
science. The NASA Scientific and Technical 
Information (STI) Program Office plays a key part 
in helping NASA maintain this important role. 

The NASA STI Program Office is operated by 
Langley Research Center, the Lead Center for 
NASA's scientific and technical information. The 
NASA STI Program Office provides access to the 
NASA STI Database, the largest collection of 
aeronautical and space science STI in the world. 
The Program Office is also NASA's institutional 
mechanism for disseminating the results of its 
research and development activities. These results 
are published by NASA in the NASA STI Report 
Series, which includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant 
phase of research that present the results of 
NASA programs and include extensive data 
or theoretical analysis. Includes compilations 
of significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA's counterpart of peer- 
reviewed formal professional papers but 

has less stringent limitations on manuscript 
length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or 
of specialized interest, e.g., quick release 
reports, working papers, and bibliographies 
that contain minimal annotation. Does not 
contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or cosponsored by 
NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific 
and technical material pertinent to NASA's 
mission. 

Specialized services that complement the STI 
Program Office's diverse offerings include 
creating custom thesauri, building customized 
data bases, organizing and publishing research 
results . . . even providing videos. 

For more information about the NASA STI 
Program Office, see the following: 

• Access the NASA STI Program Home Page 
at http:/ lzvww.sti.nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA Access 
Help Desk at 301-621-0134 

• Telephone the NASA Access Help Desk at 
301-621-0390 

• Write to: 

NASA Access Help Desk 

NASA Center for AeroSpace Information 

7121 Standard Drive 

Hanover, MD 21076 


NASA/TM— 2000-210594 


AIAA-2001-0986 



Computation of an Underexpanded 3-D 
Rectangular Jet by the CE/SE Method 

Ching Y. Loh, Ananda Himansu, and Xiao Y. Wang 
Taitech Inc., Cleveland, Ohio 

Philip C.E. Jorgenson 

Glenn Research Center, Cleveland, Ohio 


Prepared for the 

39th Aerospace Sciences Meeting and Exhibit 

sponsored by the American Institute of Aeronautics and Astronautics 
Reno, Nevada, January 8-11, 2001 


National Aeronautics and 
Space Administration 


Glenn Research Center 


December 2000 


Acknowledgments 


The authors wish to thank Dr. L.S. Hultgren for fruitful discussions. This work received support from the 
Supersonic Propulsion Technology Project Office of NASA Glenn Research Center. 


This report is a preprint of a paper intended for presentation at a conference. Because 
of changes that may be made before formal publication, this preprint is made 
available with the understanding that it will not be cited or reproduced without the 

permission of the author. 


Available from 


NASA Center for Aerospace Information 
7121 Standard Drive 
Hanover, MD 21076 
Price Code: A03 


National Technical Information Service 
5285 Port Royal Road 
Springfield, VA 22100 
Price Code: A03 


Available electronically at http: / / gltrs.grc.nasa.gov/GLTRS 


AIAA-2001-0986 


COMPUTATION OF AN UNDEREXPANDED 3-D RECTANGULAR JET 

BY THE CE/SE METHOD 

Ching. Y. Loh* Ananda Himansu,* Xiao Y. Wang; 

Taitech, Inc. 

Beaver Creek, Ohio 44135 

Philip C. E. Jorgenson 

National Aeronautics and Space Administration, Glenn Research Center 

Cleveland, Ohio 44135 


Abstract 

Recently, an unstructured three-dimensional 
space-time conservation element and solution 
element (CE/SE) Euler solver [1] was devel- 
oped. Now it is also developed for paral- 
lel computation [2] using METIS [3] for do- 
main decomposition and MPI (message pass- 
ing interface) [4]. The method is employed 
here to numerically study the near-held of a 
typical 3-D rectangular under-expanded jet. 

For the computed case — a jet with Mach num- 
ber Mj = 1.6, with a very modest grid of 1.7 
million tetrahedrons, the how features such as 
the shock-cell structures and the axis switch- 
ing, are in good qualitative agreement with ex- 
perimental results [8], [20]. 

1 Introduction 

An under-expanded supersonic jet radiates mixing 
noise, broadband shock-associated noise, as well as 
screech tones under certain conditions. These compli- 
cated and technologically important physical phenomena 
have been the topic of many experimental and theoreti- 
cal investigations, see Tam’s review papers [10, 1 1] for a 
comprehensive list of references. Generally, the mixing 
noise is directly associated with large-scale structures, 
or instability waves, in the jet shear layer, whereas the 
broadband shock-associated noise and screech tones are 
associated with the interaction of these waves with the 
shock-cell structure in the jet core. The distinct screech 
tones arise due to a feedback loop, i.e. 9 part of the acous- 
tic waves generated by the wave/shock-cell interaction 
propagate upstream and re-generate the instability waves 
at, or in the vicinity of, the nozzle lip. More details can 
be found in the review papers [9-11] and the references 
therein. 

* Member A1AA 
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Jet noise is a challenging topic in computational aeroa- 
coustics, in particular, near-field noise computation in the 
presence of shock cells in the jet core. In this situation, 
the computational scheme is required on one hand to re- 
solve the acoustic waves without introducing too much 
dispersion error and numerical dissipation, while on the 
other hand, it is required to capture shocks, or other non- 
linear phenomena, near or inside the jet correctly. In ad- 
dition, non-reflecting boundary conditions must be im- 
plemented, which is more difficult to accomplish in the 
near field than in the far field. 

The ‘Space-Time Conservation Element and Solution 
Element Method' [5-7], or the CE/SE method for short, 
is a scheme that meets the above requirements. The 
CE/SE scheme possesses attractive properties for aeroa- 
coustics computations in that: (/ ) it possesses low dis- 
persion and dissipation errors; (ii) its ‘built-in’ shock- 
capturing nature makes the computation of shock-cell 
structures simple and accurate; (iii) the non-reflecting 
boundary conditions are simple and effective and can 
be applied in the near field of the jet without introduc- 
ing excessive errors; and (/v ) the scheme accurately pre- 
dicts the vorticity, which plays an important role in the 
noise generating mechanism. A detailed description of 
the CE/SE method can be found in the reports of Chang 
et a I. [6,7]. As demonstrated in our previous papers, the 
CE/SE scheme is well suited for computing waves on 
compressible shear flows [12] as well as vorticity/shock 
interactions [12, 13], both being comer stones of the jet- 
noise phenomena. Recently, the CE/SE Euler solver has 
been developed and used for 3-D compressible flow and 
aeroacoustics computations [1]. The solver is based on 
unstructured grids, and is robust for problems of general 
geometry. In order to increase the computational effi- 
ciency and to reduce the turn-around time, the 3-D Euler 
code has been parallelized [2] using METIS [3] for do- 
main decomposition and MPI [4] for message passing 
between processors. 

In this paper, the near-field of an under-expanded 3- 
D rectangular supersonic jet is investigated numerically 
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by using a CE/SE Euler solver. The paper is arranged 
as follows: The unstructured CE/SE Euler scheme, and 
the parallel computation are briefly described in Section 
2. Section 3 illustrates the initial and boundary condi- 
tions for the 3-D rectangular jet problem, in particular, 
the novel non-reflecting boundary conditions, which are 
based on flux balance. The numerical results are pre- 
sented and compared to experimental findings [ 8 ] in Sec- 
tion 4. Conclusions are given in Section 5. 

2 The 3-D CE/SE Euler Scheme 

The 3-D CE/SE method systematically solves a set of 
integral equations derived directly from the physical con- 
servation laws, with space and time treated ‘on the same 
footing'. Because of its integral formulation, the scheme 
naturally captures shocks and other discontinuities in the 
flow. Both dependent variables and their derivatives are 
solved for simultaneously and, consequently, the flow 
vorticity can be obtained without reduction in accuracy. 
Non-reflective boundary conditions are also easily im- 
plemented because of the flux-conservation formulation. 
Details on the 3-D CE/SE Euler method can be found in 
Wang et al.[ 1]. In the following subsections, the CE/SE 
method is briefly reviewed. 

2.1 Conservation Form of the 3-D Unsteady 
Euler Equations 

Consider a dimensionless conservation form of the un- 
steady 3-D Euler equations for a perfect gas. Let p, u, 
v,w , p, and 7 be the density, velocity components in x, y 
and 2 directions, static pressure, and constant specific 
heat ratio, respectively. The Euler equations then can be 
written in the following vector form: 

U t + F x + G y + H z = 0 , ( 1 ) 

where x, y , z and t are the spatial coordinates and time, 
respectively, and the conservative flow variable vector U 
and the flux vectors, £, G and if, are given by: 
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with 

U\ = p, U 2 = pu , U s = pv, Ua = pw, 

Us = p/( 7 - 1 ) + p{u 2 + v 2 + w 2 )/2 
Fi = C/ 2 , 

F2 = (7- l)£fc+[(3 - 7)^2 - (7 - 1 )(ul + UD] /2 Uu 
F 3 = U 2 U 3 /U 1 , 

F 4 = U 2 U 4 /U 1 , 


F- 0 = 7 U 2 u 5 /u, - (7 - 1 )U 2 [Ul + 0 3 2 + Ul\ /2Ul 

Gi=U 3 , g 2 = u 2 u 3 /u 1 , 

G 3 = (7-1)^ b +[(3 - 7)08 - (7 - l)(Ui + Ul)} /2U U 
Gi = U 3 U 4 /U 1 , 

G 5 = 70305/01 - (7 - 1)03 [0 2 2 + Ul + Ul] /2 Ul 
01=04, 02 =0204/01, 

03 =0304/01, 

04 = (7-1)05+ [(3 - 7)04 2 - (7 - 1)(0 2 + 0 3 2 )] /201, 

05 = 70405/01 - (7 - 1)04 [0 2 2 + Ul + Ul] /2Ul 

By considering ( x,y,z,t ) as coordinates of a four- 
dimensional Euclidean space, £4, and using Gauss' di- 
vergence theorem, it follows that Eq. (1) is equivalent to 
the following integral conservation law: 

/ Qm • dS = 0 . m = 1 , 2 , 3, 4, 5, ( 2 ) 

JS(V) 

where S(V) denotes the surface around a volume V in 

£4 and Q m = (£ m ,G m , U m ) . 

2.2 CE/SE Structure 

In the 3-D CE/SE scheme, as in its 1 -D and 2-D coun- 
terparts, the flux conservation relation in space-time is 
the only mechanism that transfers information between 
nodes. The conservation element, C£, is the 4-D space- 
time finite volume to which the integral flux condition 
(2) is to be applied. Solution discontinuities are allowed 
to occur in the interior of a conservation element. A so- 
lution element, S£, associated with a grid node is here 
a set of seven interface hyperplanes in £4 that passes 
through this node. Note that these hyperplanes are ac- 
tually 3-D space volumes or space-time volumes. The 
solution, i.e. U , U x , U y , and U z is calculated at this 
node. Within a given solution element 5£(j, n), where 
j is the node index of the unstructured grid and n the 
number of the time step, the flow variables are not only 
considered continuous but are also approximated by lin- 
ear Taylor expansions: 

U*(x,y,z,t;j, n) = U] + (U x )^(x-Xj)+ 

(U y )?(y - Vj ) + (U 2 )^z - zj) + (U t )”(t - t n ), (3) 
F*(x,y,z,t]j, n) = F” + {F x )](x - Xj)+ 
(F y )?(y - yj ) + (F 2 )”(z - zj) + ( F t )?(t - t n ), (4) 
G*(x,y,z,t;j,n) = G] + (G x )^(x - xj)+ 
(G y )](y - yj) + (Gz)'j(z — Zj) + (G t )^{t - t n ), (5) 
H*(x,y,z,t;j,n) = if" + (H x )?(x - xj)+ 
(Hy)J(y - Vj ) + (H z )”(z - zj) + (H t )?(t - t n ), ( 6 ) 


NAS A/TM— 2000-2 1 0594 


2 

American Institute of Aeronautics and Astronautics 


where the partial derivatives of £, G and H can be re- 
lated to the corresponding ones of U by using the chain 
rule and U t can be directly obtained from ( 1 ). 

The discrete approximation of (2) is then 



Each S(CE) is made up by “surface segments” belong- 
ing to n\>o neighboring S£’s. All the unknowns are 
solved for based on these relations. No extrapolations 
(interpolations) across a stencil of cells are needed or al- 
lowed. 

The CE/SE scheme is naturally adapted to the topol- 
ogy of an unstructured grid. In the current 3-D case, as 
shown in Figure 1 , the four vertices A, B,C,D of any 
tetrahedron in an unstructured mesh form a dodecahe- 
dron along with the cell centers N1,N2, NS and 7V4 
of its four neighboring tetrahedrons. The cell center of 
the current tetrahedron ABCD is O. The dodecahedron, 
being the projection of an £4 space-time volume V in 
(2) onto the 3-D space, is formed by four hexahedrons: 
O - ABD -N1,0- BCD -N2,0- CD A - N3 
and O - ABC — N4. These hexahedrons are the projec- 
tions of space-time £4 CE’s onto 3-D space. There are 
totally 20 scalar unknowns U, U x , U y , and U z at O. 
Each of the 4 neighboring cells provides 5 scalar equa- 
tions (2), totaling to 20 equations. At each marching time 
step, the conservative flow variables at any tetrahedron 
center is updated based on the flow variables at its four 
neighboring tetrahedron cell centers of the previous time 
step. Details of the procedure can be found in Wang and 
Chang [1]. 

2.3 Running the CE/SE Euler Code on a 
Parallel Computer 

In aeroacoustic computations, due to the stringent re- 
quirement of acoustic wave resolution, a certain number 
of grid points (or cells) per wavelength must be main- 
tained. Generally, with the capability of the current com- 
puters, 2-D aeroacoustics problems can be properly han- 
dled by a single processor (CPU) without excessively 
long run times or memory problems. However, for 3-D 
aeroacoustics computations, from the viewpoint of com- 
putation turn-around time and memory sizes, parallel 
computation with multi-processors becomes necessary. 

Four steps are taken in the parallelization of the 3-D 
CE/SE Euler solver: 

1 . an unstructured grid is generated for the entire do- 
main; 

2 . the domain is decomposed into subdomains accord- 
ing to the assigned number of processors (CPU’s), 
using the METIS code. METIS is a mesh partition- 
ing code in the public domain and available from the 


University of Minnesota [3]. It can also be freely 
downloaded from the web. 

3. the flow is computed in each subdomain with the 
corresponding CPU using the CE/SE Euler solver. 
The CE/SE code is modified using the MPI library 
calls. CPU assigned to one subdomain may commu- 
nicate with other CPU’s assigned to the neighboring 
subdomains and exchange their computed results. 
MPI (message passing interface) is an interproces- 
sor communication protocol standard. Public do- 
main implementations of MPI as a software library 
package are prepared by the Argonne National Lab- 
oratory. 

4. at the end of computation, the subdomain data is 
recombined. 

Figure 5 is an illustrative description of the above 
steps. Detailed implementation of the parallelization of 
the CE/SE method can be found in Himansu et a/.[2]. 

3 The 3-D Rectangular Jet Problem 

Consider a rectangular jet as sketched in Figure 2. The 
aspect ratio of the nozzle is 5. The nozzle lip is allowed 
to have a finite thickness (3 cells) to complete the feed- 
back loop for the possible self-sustaining screech tone. 
The jet Mach number, Mj, is 1.6, representing an under- 
expanded off-design status. The ambient flow around the 
jet is stationary. These conditions correspond to the ex- 
perimental setup of Raman [ 8 ]. 

In this investigation, our attention is focussed on the 
near field of the jet nozzle since this is the source re- 
gion for aeroacoustic noise. The inner short side width, 
£, of the jet nozzle is chosen as the length scale. The 
computational domain is a rectangular hexahedron with 
length, L, of 1 6 £, width, W, 14 £, and height, £, 5.6£ 
as shown in Figure 2. The nozzle is located in the central 
area of the computational domain and protrudes into the 
domain with l =2D (Figure 2). The thickness of nozzle 
lip wall is 2 cells in order that a screech feedback loop 
can be formed. The unstructured tetrahedron grid cur- 
rently used is generated by cutting a hexahedron cell into 
six tetrahedrons as shown in Figure 3. Currently, includ- 
ing the nozzle interior, there are 80 x 60 x 60 such hex- 
ahrdron cells in the computational domain, which form 
about 1.7 million tetrahedrons. These hexahedron cells 
are non-uniform in size in order to have better resolution 
around the jet shear layers. The computational scheme is 
of the a-e type [6,7] with a = 0 and e = 0.5. 

3.1 Initial Conditions 

As shown in Figure 2, the convergent nozzle protrudes 
into the computational domain so there is some room to 
capture any waves that propagate upstream from the noz- 
zle exit. A conceivable initial condition is to assume that 
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the pressure is uniform in the entire computational do- 
main and equal to the ambient pressure, except at the jet 
nozzle exit where a higher pressure p e is imposed as the 
boundary condition since the jet is under-expanded. We 
will also make the reasonable assumption that the tem- 
perature in the plenum equals the ambient one for this 
cold jet. 

Let the density, axial velocity, and temperature in the 
jet core, be used to scale density, velocity and tempera- 
ture respectively. Hence, the initial conditions in the jet 
core are given by 


Pj = 1» Tj = 1, Uj = 1, vj = 0, Pj = 


7 MV 


where T denotes the temperature. The initial conditions 
in the ambient region are given by 


Pa = 1 /T a , T a = 1 + |( 7 - l)Mj, 


Type I and Type II CE/SE non-reflecting boundary con- 
ditions as described in the next subsection are imposed 
respectively. 

3.3.2 Non-Reflecting Boundary Conditions 

In the CE/SE scheme, non-reflecting boundary condi- 
tions (NRBC) are constructed so as to allow fluxes from 
the interior domain of a boundary CE to smoothly exit 
the domain [18]. There are various versions of the non- 
reflecting boundary condition, and in general they have 
proven to be well suited for aeroacoustic problems [ 1 2- 
17]. The following are the ones employed for the 3-D 
problem: 

1. Type I: For a ghost grid node (— j, n) lying outside 
the domain boundary and being a mirror image of 
an interior node (j, n), where j is the node index 
number and n the number of the time step, the non- 
reflecting boundary condition (type I) requires that 


u a = 0, 


v a = 0, 


Pa ~ 1 M] ’ 


where the result for the temperature follows from the as- 
sumptions of constant total enthalpy in the jet flow, and 
the plenum and ambient temperatures are equal. 

By using the assumption of constant total enthalpy and 
the condition for isentropic flow, it follows that the pres- 
sure p e at the nozzle exit can be related to the pressure 
Pj in the jet core by 


p 1= [ 1 + |(7-1 )M] Y^ 

Pj [l + i( 7 -l )Mf_ 

Here, M e = 1 is the Mach number at nozzle exit when 
the flow is choked. 

3.2 The jet shear Layer 

In the present numerical computations, no jet shear layer 
exists initially in the field. Since the pressure at the noz- 
zle exit is always higher than the ambient one, the flow 
pushes its way in the ambient atmosphere and forms the 
jet and jet shear layers. The spreading of the jet shear 
layer is caused by the mixing and momentum exchange 
in the shear layer and is very important for the noise gen- 
eration (e.g. [11]). 

3.3 Boundary Conditions 

3.3.1 Inflow Boundary Condition At the in- 
let boundary, the dependent variables and their spatial 
derivatives are specified to be those of the ambient flow, 
while at the nozzle exit, the elevated pressure p = p e is 
imposed, i.e. the jet is under-expanded, as in the physical 
experiment. 

No artificial forcing of the shear layer is imposed at the 
nozzle. At the surrounding and outflow boundaries, the 


(U X ) n _J = (Uy) n _J = ( U z ) n _J = 0, 

while U n _j is kept fixed at the initially given steady 
boundary value. 

2. Type II: At the downstream boundary, where there 
are substantial gradients in the y , z directions, the 
non-reflective boundary condition (type II) requires 
that 

(U x )!j = 0 , 

while U^j, (U y )*j and (U z )”_j are now defined 
by simple extrapolation from the interior, i.e., 

u n -j = u?-\ 

(Uy) n _J = (t/,)?- 1 , ([/;)-, = (IMr 1 - 

These non-reflecting boundary conditions (NRBC) are 
consistent with their counterparts in lower dimensions 
(1-D and 2-D). They are simple to implement and effec- 
tive as well. 


4 Numerical Results 

In this section, the numerical results for the 3-D under- 
expanded rectangular jet are presented and compared to 
experimental results [8]. At the moderate supersonic jet 
Mach number Mj = 1.6, the overall motion in the ex- 
periment [8] is in truly 3-D mode. For imperfectly ex- 
panded jets such as this one, a quasi-periodic shock-cell 
structure is formed in the jet plume. The streamwise- 
growing instability waves naturally occurring in the jet 
shear layer interact with this shock-cell structure and 
generate broadband shock-associated noise and under 
certain conditions screech tones develop through a feed- 
back loop. 
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It is important to realize that even if no harmonic forc- 
ing is intentionally imposed on or introduced in the nu- 
merical simulation, the jet shear layer is in actuality con- 
tinuously stimulated at a very low level as a result of trun- 
cation, round-off, and discretization errors (all of which 
can be characterized as numerical noise — and, in a sense, 
are analogous to environmental background noise in ex- 
periments). These growing perturbations interact with 
the shock-cell structure of the jet plume and thereby gen- 
erate acoustic waves. 

4.1 Test Run of Parallel Computation 

Before commencing parallel computation of the rectan- 
gular jet problem, the numerical procedure is first tested 
on a simple jet problem with fewer grid nodes, so that it 
can also be run serially on a single CPU with reasonable 
turn-around time. 

Figure 4 depicts a sketch of the simplified jet problem. 
The computational domain is a rectangular hexahedron. 
At the inlet plane, boundary conditions are set to the val- 
ues of the stationary ambient flow except at the center 
square jet nozzle, a jet flow Mj — 1.2 is imposed. The 
stationary ambient flow is also used as the initial condi- 
tions. Type I and II NRBC’s are applied accordingly. A 
modest grid of 30 x 30 x 30 rectangular hexhedron blocks 
is used to generate the unstructured grid. The problem 
is then run first in parallel mode on an SGI workstation 
cluster using 8 processors of 195 MHz, and second, seri- 
ally on a PC with Pentium III 600MHz CPU for 20,000 
time steps. In both cases, only single precision arith- 
metic is used. Figure 5 illustrates the schematic diagram 
of the parallel computation. The computed results are 
practically identical, as demonstrated in Figure 6. While 
the computation on the serial machine takes about 9 sec- 
onds CPU time per time step, it only takes about 1 second 
(CPU time) per time step in the parallel mode on the SGI 
workstation cluster. The advantage of parallel computa- 
tion is thus clearly demonstrated. 

4.2 The Rectangular Jet Problem 

The grid, initial and boundary conditions for the cur- 
rent rectangular jet problem are already described in the 
previous section. The actual computation was run for 
60,000 time steps in order to ensure that all transients are 
convected out of the computational domain. 

In the preprocessing stage, domain decomposition 
with METIS and the data setup for each processor only 
takes a few minutes to complete. Sixteen processors on 
an SGI workstation cluster are employed for the com- 
putation. These CPU’s are R 10000 chips ranging from 
195MHz to 250 MHz in speed. With these 16 processors, 
the average CPU time consumption is about 5.1 second 
per time step. However, depending on the load of the SGI 
cluster, it takes about 5-6 days turn-around time to run 
60,000 time steps, which is still acceptable in our view. 


In the following, snapshots of flow variables are plot- 
ted and compared to available experimental results or ob- 
servations. It should be emphasized that no buffer zone 
at the domain boundaries is used in these computations. 
The CE/SE NRBC is robust enough to yield a solution 
free of noticeable spurious reflections. 

4.3 Jet Shock-Cell Structure 

Experimental results for supersonic jets are often doc- 
umented in terms of Schlieren pictures. It is straight- 
forward to construct Schlieren graphs from the numeri- 
cal results without loss of accuracy in view of the nature 
of the CE/SE method, since density derivatives are di- 
rectly solved for as unknowns. Figure 7 shows a compar- 
ison of numerical Schlieren (density gradient modulus) 
contours from the current computations and the experi- 
mental Schlieren photograph in [8]. The numerical result 
is a snapshot of the instantaneous density gradient con- 
tour plot on the mid sectional plane of the narrow side 
at the end of the computation (60,000 time steps). The 
shock cells that have developed so far in the simulation 
agree quite well with the experimental ones. In the same 
Figure 7, numerical Schlieren on the mid sectional plane 
of the wide side is also shown. Shock-cell structures in 
the jet plume are clearly demonstrated. 

Similar shock-cell structures are also observed from 
the pressure, u- velocity, and Mach number contours in 
Figures 8, 10, and 11. 

4.4 Near-Field Jet Flow 

In the computation, no artificial forcing is imposed at the 
nozzle exit at all, yet the shock-cell structure and insta- 
bility waves seem to be sustainable (at least to the stage 
of 60,000 time steps). This implies a feedback loop is 
likely formed around the shear layer between the noz- 
zle lip wall and the shock-cell structure, although more 
refined grid and much more time steps are required to 
verify this. (Based on our experiences with 2-D under- 
expanded axisymmetric sonic jet computation, screech 
wave takes place when the grid is fine enough. ) 

Figure 8 shows the isobars on the mid sectional planes 
of both the narrow and the wide sides. Similar shock- 
cell structures as in Figure 7 are found. Pressure iso- 
surfaces at two close but different levels are plotted in 
Figure 9. It is clearly demonstrated that these pressure 
waves at the two edges of the wide side jet plume form a 
helical pattern. Qualitatively, this pattern agrees with the 
observation of Quinn [20] and with the observation of 
Westley and Woolley [19] on screech of spinning mode. 
Although these waves look like Mach radiation, at the 
present stage, with a very coarse grid of about 1.7 mil- 
lion cells, it is not sure if these pressure waves are truly 
radiating Mach waves. Investigation with a finer grid is 
needed. 
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Figure 1 0 illustrates the u velocity contours at the mid- 
sectional planes of the narrow and wide sides. The out- 
skirts of the jet plume on the narrow side form barrel-like 
shapes according to the shock-cell structure and spread 
out due to mixing. This is in contrast to the axisymmet- 
ric mode, where the outskirts of the jet plume form only 
smooth lines. Another interesting observation from the 
numerical results is the “axis switching" phenomenon, 
the jet core seems to keep reducing its width on the wide 
side and increasing its width on the narrow side due to 
shear layer spreading. Therefore, a few diameters ( D ) 
downstream from the nozzle exit, the jet core appears to 
be roundish, although at the nozzle exit it is rectangular 
with an aspect ratio of about 5. The same observation 
can also be found in Figure 1 1 . 

5 Concluding Remarks 

In this paper, the recent unstructured CE/SE Euler 
scheme is extended to 3-D computation for a super- 
sonic under-expanded rectangular jet. The state-of-the- 
art techniques of METIS and MPI are applied in the par- 
allelization of the unstructured CE/SE Euler code, and 
proved to be effective and efficient. 

The advantages of the CE/SE scheme found previ- 
ously in 2-D computations are confirmed again for 3-D 
simulations: 

1 . The method is robust and the implementation is ‘ef- 
fortless’ in that no special treatment and parameter 
selections are needed; 

2. The NRBC is simple and effective, no buffer zone 
is used; 

3. The method handles both linear and nonlinear prob- 
lems and is particularly advantageous for near field 
jet-noise computation; 

4. Unstructured grid and parallelization help to handle 
large scale problems with complicated geometry. 

With a very modest unstructured grid of about 1.7 mil- 
lion elements (cells) and single precision computation, 
many aspects of the computed results are in good qualita- 
tive agreement with experimental findings [8], [20]. With 
1 6 processors in an SGI workstaion cluster, the CPU time 
consumed per time step is about 4 seconds. However, a 
self-sustained oscillation, i.e. the screech tone, is yet to 
be achieved in our simulation. At the nozzle lip, a finite 
lip wall is already added to facilitate the forming of the 
feedback loop. Based on current encouraging results of 
2-D axisymmetric jet screech tone computation, we are 
sure the screech can be computed with refined grid. Also, 
in order to achieve a more physically relevant simula- 
tion, a Navier-Stokes CE/SE code with LES (large eddy 
simulation) or turbulence modeling capability is need to 


account for the strong momentum exchange in the shear 
layer. The results will be reported in the future. 
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D 



A, B, C, D - tetrahedron vertices, 
O cell center, N1,N2,N3,N4 - 
neighboring cell centers 



showing part of the hexahedron 
O-ABD-Nl (shaded area) 

Figure 1: Tetrahedron and dodecahedron in an unstruc- 
tured grid. 
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Figure 4: Sketch of the test problem. 



Figure 2: Sketch of the rectangular jet problem. 
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Figure 3: Partition of a hexahedron into 6 tetrahedrons. 


Figure 5: Schematic diagram of parallel computation 
(showing 5 processors). 
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Figure 6: Comparison of numerical results by parallel and serial computations. 
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Figure 7: Comparison of numerical Schlieren graph (instantaneous snapshot) and the experimental time averaged 
Schlieren graph ; showing good agreement of the first few shock cells. 
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isobars, narrow side mid sect, plane 



isobars, wide side mid sect, plane 

Figure 8: Pressure contours on the narrow and wide side 
mid sectional planes. 



Figure 9: Pressure iso-surfaces at 2 different levels, 
showing the shock-cells, the helical Mach radiations. 



u-contours, wide side mid sectional plane 


Figure 10: u- velocity contours 



M, narrow side mid sectional plane 



Figure 1 1 : Mach number contours. 
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